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ABSTRACT 

We present an algorithm for generating merger histories of dark matter haloes. The 
algorithm is based on the excursion set approach with moving barriers whose shape 
is motivated by the ellipsoidal collapse model of halo formation. In contrast to most 
other merger-tree algorithms, ours takes discrete steps in mass rather than time. This 
allows us to quantify effects which arise from the fact that outputs from numerical 
simulations are usually in discrete time bins. In addition, it suggests a natural set of 
scaling variables for describing the abundance of halo progenitors; this scaling is not as 
general as that associated with a spherical collapse. We test our algorithm by compar- 
ing its predictions with measurements in numerical simulations. The progenitor mass 
fractions and mass functions are in good agreement, as is the predicted scaling law. 
We also test the formation-redshift distribution, the mass distribution at formation, 
and the redshift distribution of the most recent major merger; all are in reasonable 
agreement with N-body simulation data, over a broad range of masses and redshifts. 
Finally, we study the effects of sampling in discrete time snapshots. In all cases, the 
improvement over algorithms based on the spherical collapse assumption is significant. 
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• 1 INTRODUCTION 



Most models of galaxy formation in a hierarchical universe 
assume that the merger history of the surrounding dark mat- 
ter halo plays an im portant role in determ i ning the proper- 
ties of a galaxy (e.g. IWhite fc Ree3 (|l978l ) ; IWhite fc Freii^ 
(|l99ll '): lBaughl l|2006l 'l and references therein). Although halo 
merger histories can be measured using N-body simulations, 
these can be time co nsuming and computationally intensive 
l|Springel et al.ll2005h . This has fueled considerable study of 
the formation and merger histories of dark matter haloes 
from a Monte Carlo perspective. Monte Carlo merger trees 
have the advantage of being fast and one may easily probe 
mass regimes inaccessible to current N-body simulations. 
Moreover, unlike N-body experiments, the cosmology and 
initial conditions may be easily modified^ 

The excursion set framework (IBond et al.l Il99ll : 
1 1 — — ' ^ 

acev fc Cole 1993), which is motivated by the pioneering 

work of lPress fc Schechter ( 1974 ) , provides the basis for cur- 
rent models of halo assembly. Initially, this framework was 
based on the assumption that haloes form from a spher- 
ical c ollapse, of the type first described by ICunn fc GottI 
Fast algorithms for generating halo merger trees, in 



which haloes were assumed to form from a spherical cqL 
lapse , were developed in the 1990s ( KauflFmann fc Whitt 
19931: ISomerviUe fc Kolattl Il999l: jSheth fc LemsonI 1199^ 
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Cole et al.l l2000l ) (see IZhang et al.l (|2008bh for a review) . 
However, spherical collapse overpredicts (underpredicts) 
the abundance of haloe s in the low (high) mass regime. 
To address these issues, ISheth fc To"rm on (1999) extended 
the excursion set framework to include ellipsoidal c ollapse 
ISheth. Mo fc TormenI [2OO1I : ISheth fc Tormenl l2002l 'l. This 
clearly showed that merger-trees which assume spherical col- 
lapse are inadequate. 

iHiotelis fc Del Popold (|2006l ) describe a merger tree al- 
gorithm which extends some of the older algorithms to incor- 
porate aspects of the ellipsoidal collapse results. In addition, 
a number of new algorithms have recently b een published 
ijParkinson et al.ll2008l : lNeistein fc Dekei2008h : although ef- 
ficient and accurate, such methods side-step the idea of el- 
lipsoidal collapse altogether. Moreover, these methods are 
callibrated to match N-body simulations, and are therefore 
limited by the accuracy and scope of these simulations. 

The aim of the present work is to provide a merger 
history tree algorithm which is based explicitly on the 
excursion-set formalism with ellipsoidal collapse. The most 
significant difference between the algorithm we derive and 
all the others described above is that it takes discrete steps 
in mass rather than time. This feature allows us to study a 
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few problems which are more difficult to address with the 
other methods. 

After we completed this project, IZhang et al] (|2008bl ) 
presented an alternative algorithm with ellipsoidal collapse 
(and discrete-time snapshots). This algorithm generates pro- 
genitors across many mass bins and then assigns them 
to final haloes. In t his s ense, it is very similar to the 
iKauffmann fc White! (|l993l ') merger tree, but solves many of 
its problems. Although this tree is quite difi'erent from ours, 
it also attempts to describe the merger history of a halo from 
the excursion-set formalism. These two approaches show 
that generating merger trees without tuning to N-body sim- 
ulations is a quite non-trivial, yet interesting, challenge. 

A review of background material and a description of 
our algorithm are given in Section [J] Section |3] compares our 
algorithm with excursion-set theory predictions, and with 
measurements in N-body simulations. The tests include the 
progenitor mass fractions and mass functions, the formation- 
redshift distribution, the mass distribution at formation, and 
the redshift distribution of the most recent major merger. 
Section [4] summarises our findings and discusses possible ap- 
plications of our algorithm, an outline of which is in Ap- 
pendix [X] 



2 MERGER TREES IN THE EXCURSION SET 
APPROACH 

In the excursion set approach, the problem of estimating 
halo abundances is mapped to one of estimating the distri- 
bution of the number of steps a Brownian-motion random 
walk m ust take before it first crosses a barrier of specified 
height ijBond et al.lll99ll l. In this approach, the height of 
the barrier plays a crucial role. 



2.1 Constant and moving barriers 

The Press-Schechter mass function is associated with barri- 
ers of constant height - such barriers arise naturally in mod- 
els in which haloes form from a spherical collapse model. In 
constrast, in the ellipsoidal collapse model, barriers of the 
form 



S 



(1) 



are more natural. Here Sc is the overdensity required for 
spherical collapse - it is a monotonically increasing function 
of redshift, and it is given by 5co/D{z), where Sco = Sc{z — 
0) ~ 1.686 and D{z) is the linear growth factor. S is a 
monotonically decreasing function of halo mass, given by: 



(2) 



S{m) = a(m) = — / dfc fc^ P(k) W^(kR) 



where P{k) is the initial power spectrum of density fluc- 
tuations, linearly evolved to the present time, W is the 
Fourier transform of W{x) = (3/a::"^)(sina; — a; cos a;) , R = 
(3rn/47rp)^^^, and p is the comoving background density. At 
large R, the overdensity contained in the associated volume 
is practically zero. As R (and m{R)) decreases, S{R) in- 
cr eases, and 5r execute s a random walk. We refer the reader 
to lLacev fc Cold (|l993l ) for more details. 

Consider a barrier B{S,5c{z)), as in equation ([T]), and 



an ensemble of random walks which start from the origin: 
{S,5) = (0,0). The excursion set approach maps the distri- 
bution of S, the number of steps a random walk must take to 
first cross such a barrier, to the fraction of mass in m-haloes 
at redshift z. This quantity is associated with the so-called 
unconditional mass function. The conditional mass function 
of high redshift progenitors of a more massive final A/-halo 
at some lower redshift Z is modeled using walks which start 
from {Sm , 5c{Z)) instead. 

The shape of the mass functions (unconditional and 
conditional) is determined by the shape of the barrier, en- 
coded in the (q, /3,7) parameters. The spherical collapse 
model is associated with ((?,/3, 7) — (1,0,0), whereas ellip- 
soidal collapse has (0.707,0.47,0.615). When 7 > 1/2, not 
all w alks are guaranteed to c ross the ellipsoidal collapse bar- 
rier (ISheth fc Tormenll2002l ). Moreover, the barriers associ- 
ated with two different times may intersect; of course, this 
never happens for the spherical collapse barriers. Sheth & 
Tormen suggest that this intersection of barriers may repre- 
sent the possibility that haloes can fragment. This compli- 
cates the algorithm we describe below, so we instead study 
the limiting case of 'square-root' barriers for which 7 — 1/2: 



(3) 



The predicted halo abundances associated with {q, j3, 7) = 
0.55, 0.5, 0.5) are very s i milar to those in simulations 
Mahmood fc Raieshllioosi : iMoreno et. al.ll2008h . Moreover, 

the first cross ing distribution of this barrier is known 

(jBreimanlll9"67h . 



2.2 A conditional scaling symmetry 

Recall that the unconditional mass function is associated to 
the first crossing distribution associated with walks which 
start from the origin: (5*, 5) = (0, 0). When constant barriers 
are use, it can be expressed self-similarly as 



f{m\z)dm = f{S\5c)dS = f{u)du, 



(4) 



where u = S'^/S. The conditional mass function of m-haloes 
at z that end up in bound objects of mass M > m a,t Z < z 
is given by f{m,z\M, Z)dm = f{Sm, Si\Sm , 5o)dSm, where 
where Si = 5c{z), So = 5c{Z). In other words, the conditional 
mass function is associated with the first crossing distribu- 
tion of a barrier of height Si by random walks with origin 
at {Sm,So)- Because a straight-line is straight whatever the 
origin of the coordinate system, the conditional mass func- 
tion, in the spherical collapse model has the same functional 
form as that of the unconditional mass function, provided 
one sets v — {Si — So)^ /{Sm — Sm)- 

However, for the square-root barrier, a walk which 
starts from (y^^o + PV Sm, Sm) must cross a barrier of 
shape ^{Si — So) + P\/Sm — Sm + Sm- This is not quite 
of the same form as equation Q. As a result, the condi- 
tional mass function is not simply a rescaled version of the 
unconditional one. Rather, in this model. 



/(m, z\M, zo) dm = f{S,n/SM\v) d{Sm/SM), (5) 



where 



<5i — So 



(6) 



(see iBreimanI ([1963) or ICiocoh et all |20o3) for the exact 
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Figure 1. A random walk and its associated mass history. The 
dark filled circles represent the history of a halo of mass M at 
redshift z = 0. A merger {m',m — m') — > m at redshift z is 
depicted by the Sm — * S^i jump at height 5c{z). A new branch 
associated with (m — m') is connected at {S'„„,„/ , <5c(z)). The 
light-filled circles denote the mass history of this object. 



expressions). Thus, final halos of different masses will have 
similar progenitor mass functions when expressed in terms of 
Sm/SM, provided they have similar values of ry. While this 
scaling is like that for the constant barrier model, in the 
square-root barrier model, the progenitor mass function is 
not a function of the combi nation = rj^ /(Sml Sm — !)• 
This is interesting, because ISheth fc TormenI (|2002l ) have 
shown that the conditional mass function in simulations 
is not well-fit by a function of v. In what follows, we will 
present evidence that it is, however, a function of r) and 
Sm/SM separately, so the qualitatively different scaling as- 
sociated with the square-root barrier is indeed seen in sim- 
ulations. 



2.3 Mass histories and merger trees 

Figure [T] illustrates how the mass growth history of an ob- 
ject is encoded in the excursion set approach if objects form 
from a spherical collapse (also see Figure 1 of lLacev fc Cold 
l|l993l ll. The jagged line shows a random walk which starts 
from the origin: (5, 5) — (0, 0). Imagine drawing a horizontal 
line with height 5co ~ 1.686 and marking the smallest value 
of S at which the walk intersects this 'barrier' of constant 
height (5cO corresponds to the present time and Sc > Sco 
corresponds to higher redshifts). The dotted horizontal line 
denotes such barrier. Then increase the height of this bar- 
rier, and record how this value of 5* changes as S increases. 
Such mass history points are depicted as dark filled circles 
in Figure [T] The dashed lines show that S will occasionally 
jump from a small value to a larger one. Since 5* is a proxy 
for mass, and 5c for time, such a jump is a proxy for an 
instantaneous change in mass: a merger. Note that the ran- 
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Figure 2. The same random walk as in Figure [T] but now with 
square-root rather than constant barriers, illustrating that the 
mass accretion history depends on the barrier shape. In our algo- 
rithm, the new object with mass (m — m') is now connected at 



dom walk steps under such jumps are not part of the mass 
history (e.g., the gray portion with Sm < S < Sm')- 

The key to our merger tree algorithm, which is described 
in detail in Appendix El is to recognize that these jumps 
mean that there are a set of other walks which one might 
associate with this one - one for each jump. One such walk is 
illustrated by the second jagged curve, which starts at about 
the middle of the panel. If the jump from Sm to Sm' occurred 
when the barrier height was Sc{z), then this other walk starts 
from {Sm-m' ,Sc{z)). The 'merger history' associated with 
this new branch is represented by the light-shade filled circles 
in Figure [T] For every such jump, a new random walk must 
be drawn. For each jump within each of those new walks, 
the same process applies - more walks must be drawn. In 
summary, the bundle of such walks encodes the entire merger 
history of a present-day object. Notice that jumps can occur 
at any z - there is no constraint that they happen at discrete 
times. However, if one is interested in the mass function of 
progenitors at some fixed z, one simply reads-off the list of 
values of 5* at which this bundle of walks first cross Sc{z). 

So far, we have discussed how to generate trees in 
the spherical collapse model. Figure [2] shows the same 
walk as before, but now the mass growth history associ- 
ated with the walk is given by its intersection with square- 
root barriers of gradually increasing height. This shows 
clearly that the jumps in mass, and the times at which 
they occur, are modified. But the overall logic remains the 
same. Each jump gives rise to a new walk that starts from 
{Sm-m' , B{Sm-m' , Sc{z))), where B is given by equation Q. 

The natural generalisation to spherical collapse is 
to incorporate the original 7 > 1/2 barrier of 
ISheth. Mo fc TormenI (HflOl). Such a choice would compli- 
cate this algorithm significantly. First of all, as Figure [2] 
illustrates, the shape of the square-root barrier remains the 
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Figure 3. The progenitor mass fraction (left) and mass function (rigiit) at rcdsiiifts z = (0.5, 1, 2, 3, 5), for lialocs of mass M/M* = 0.06 
at 2 = 0. Filled circles show measurements in the GIF2 simulation, and open circles are from our square-root trees. The smooth solid and 
dashed curves show the exact square-root barrier solution, and the series approximation, respectively. The long-dashed curve shows the 
ellipsoidal collapse model with 7 > 1/2, and the short-dashed curve is the constant barrier prediction. Values of the scaling parameter -q 
(equation (6)1 are also shown (see Figure |6}. 



same with difFerent redshifts, except for an overall vertical 
shift. This is not the case with 7 > 1/2. As 5c increases 
(increasing redshift), the term S\^^~'^~' makes the barrier in 
equation ([T| increase less rapidly with S. A consequence of 
this is that the barriers associated with different redshifts 
cross. In the absence of crossing-barriers (e.g., constant and 
square-root barriers), one may uniquely map any point in 
the {S,5) plane to {m,z). The crossing of barriers invali- 
dates this property, making the identification of jumps with 
mergers at a given time ill-defined. 



3 COMPARISON WITH SIMULATIONS 

In this section, we compare the statistical properties of 
our merger history trees with expectations from the ex- 
cursion set theory which they are supposed to reproduce, 
and with measure ments in the GIF4J N-body simulation 
(iGao et al ] |2004bh . The simulation followed the evolution 
of 400^ particles in a periodic cubic box 110/i~^Mpc on a 
side in a fiat ACDM background cosmology with parameters 
(nm,(T8,h,nth^,n) = (0.3,0.9,0.7,0.0196,1). Fifty simula- 
tion snapshots were output, equally spaced in log(l + z). 



German Israel Fund. 
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Figure 4. Same as Figure [31 but with M/M* = 0.6. 



At each snapshot, haloes were identified using the spherical 
ov erdensity c r iterio n, adopting for virial mass the definition 
of lEke et aD l)l996l l (i.e., with virial density at ~ 324p at 
redshift zero). The particle mass is nip — 1.73 X lO^ft-^M© 
and only objects with at least ten particles are considered. 
M~k{z), defined by S'^(z) = S'(M*(z)), is the typical mass 
scale at redshift z. It is common practice to express halo 
masses in terms of M* — A'h{z — 0) (the z-dependence 
is suppressed for the present time). For this cosmology 
and initial power spectrum, = 8.7 x 10^^ h~^MQ ~ 
SOSOrrip. The simulation data and halo catalogues are 
avai lable at http://www.mpa-garching.mpg.de/Virgo See 
iGiocoli et al . (2008a) for more details regarding the post- 
processing of the simulation. 

To compare our merger histories with those in the GIF2 
simulation, we generated 2000 realisations of our tree for 
each final halo mass bin AI of interest. In all cases, the 



minimum mass considered was mdust = Af/1000, and the 
merger histories of haloes with mass below this were not fol- 
lowed (we call this minimum mass the 'branching-mass res- 
olution'). We used random walks with 10^ steps in between 
Sm and Sdust to ensure that the mass change between the 
steps was less that mdust. Having a small step size is es- 
sential to faithfully reproduce random walks in a computer. 
Moreover, if the step size is too large, we run the risk of 
missing branches. Numerical tests showed that the outputs 
converged to our results for small step sizes. See the Ap- 
pendix for more details on the implementation of our Monte 
Carlo tree. Recall that our tree does not take discrete steps 
in time. Nevertheless, for fair comparison with the measure- 
ments from the GIF2 simulation, the tree data were stored in 
the same discrete redshift bins as were output from the sim- 
ulation. We use 'Cont' to denote the original tree data and 
'Snap' for the data stored in redshift snapshots. Sections l3.3l 
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m/M m/M 

Figure 5. Same as Figures |3] but with M/M* = 6. 



and 13.41 study some merger-related quantities which are sen- 
sitive to the differences between these two ways of storing 
trees (Figures El [9] and [10}. 

3.1 The progenitor mass function 

Figures |3]l5] show the progenitor mass fractions and mass 
functions at five different redshifts {z — 0.5,1,2,3,5), for 
haloes identified ai z — with final masses given by 
M/Mi, — 0.06, 0.6 and 6. The corresponding values of rj 
(equation [Sjl are shown in each panel. In all three figures, 
filled circles show measurements in the GIF2 simulation, 
and open circles show results from our square-root trees. 
We probe the m < nip regime with our trees to verify con- 
sistency with analytic excursion-set predictions (the smooth 
curves in all the panels). The expression s we u se are given 
explicitly in Appendix A of lGiocoli et al.l l|2007l ) . The short- 



dashed curve shows the constant barrier ( /?, q) = (0, 1) 
prediction associated with spherical collapse (iLacev fc Cold 
1993). The solid curves show the exact square-root bar- 
rier solution with (/3, 5,7) = (0.5,0.55,0.5); this is a com- 
plicate d affair, i nvolv ing sums of Parabolic Cylinder func- 
tions (jBreimanl 1 19671 ). Dashed curves show the consider- 
ably simpler ap proximation to the solution which is due to 
ISheth fc Tormen (,2002 ) ; this approximation is excellent over 
the entire range of interest. The long-dashed curves show 
this same approximation for the ellipsoidal collapse barrier: 
(/3,q,7) = (0.707,0.47,0.615). The square-root barrier pre- 
diction agrees well with the 7 — 0.615 curve, except in the 
high-mass regime. This discrepancy becomes evident when 
77 > 1 and it is amplified with increasing rj (see below). 

Before we ask how our merger tree algorithm compares 
with simulations, we note that it produces progenitor mass 
functions that are well-described by the theory curves over 
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Figure 6. Th e ?7-symmetry. Dif f erent combinations of M and z with similar rj (equation [6] and Table |3!T|| . N-body simulation measure- 
ments and the ISheth fc TormenI 1I2OO2I') result with 7 > 1/2 are shown. 



a wide range of masses and redshifts. At high redshifts, 
our tree data he shghtly below the theory curves at both 
high and low m/M, and slightly above in between, although 
where the cross-over points occur depends on z and Af . In all 
other regimes, our Monte Carlo trees match the square-root 
barrier predictions. Any additional disagreement with the 
GIF 2 simulation measurements (compare open and filled 
circles) is due to limitations of the 7 = 1/2 model. 

Finally, recall that the square-root and constant bar- 
rier models make specific predictions for how the conditional 
mass functions should scale with final halo mass and time. 
Table [ST] lists pairs with similar ry, yet quite distinct values 
of M and z. Figure [S] compares the associated conditional 
mass functions. The black squares and long-dashed lines re- 
fer to the left-hand side of Table 13.11 (low final masses), 
whereas the gray triangles and short-dashed lines refer to 
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— 6 
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Table 1 . The r;-symmetry (equation used for the comparisons 
shown in Figure [6] 



the right-hand side (high final masses). Notice that the 
curves are remarkably similar to one another, as are the 
symbols. This is true despite the fact that the values of 
rj are not perfectly identical, and that f{m,z\M,Z)dm = 
f{S,n/SM\v)<i{Sr„/SM) ~ f{m/M\r,)d{m/M). The resuhs 
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Figure 7. Distribution of formation rcdshifts. Filled circles show simulation data, open circles and triangles show results from the 
square-root and constant barrier trees. Smooth curves show equation (0 with g = 1, 0.707 and 0.55 (short-dashed, long-dashed, and 
dashed), corresponding to the predicted distribution for constant barriers (spherical collapse), moving (ellipsoidal collapse) and square- 
root barriers. 



for low-mass haloes (black squares) are truncated at higher 
m/M than they are for larger M (gray triangles), simply 
because only haloes with at least ten particles are consid- 
ered. Evidently, the conditional mass functions are indeed 
functions of 77 and Sm/SM separately, rather than of the 
combination v. 



3.2 The distribution of formation redshifts 

Following iLacev fc Coll (|l993h . a halo is said to have 
'formed' when it first acquires half of its final mass. For 
a given halo mass, there is a distribution of formation red- 
shifts. This distribution is expected to peak at earlier times 
for lower mass haloes. The excursion set model with constant 



barriers provides a simple expression for this distribution of 
formation times: 

p(ljf) dtJF = 2ajF erfc(cjF/\/2) Aluf, (7) 

where 

_ Sc{zy) - 5c{za) ?7 

lof = — , = ~ — , , (ol 

y/S{M/2) - S{M) ^S{M/2)/S{M) - 1 

with 77 given by equation ((SJ. 

The filled circles in Figure[7]show the formation redshift 
distributions for haloes with masses in the range 0.9M 5C 
M ^ I.IM with logio(M/M.) = {1, . . . , -1.5} in steps of 
—0.5 in the GIF2 simulation. We use the first snapshot when 
at least half the mass is in a single progenitor as the forma- 
tion time, and make no attempt to interpolate our simu- 
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Figure 8. Distribution of the mass at formation for several final masses. The left half of each panel shows the mass just prior to 
formation, whereas the right half shows the mass just after formation. Filled circles show simulation data, open circles and triangles 
are from the square-root and constant barrier trees. The solid curve shows /ig(^) (right half) and fJ.p{fi) (left half) (equations 1101 and l9l 
respectively). The dashed curve shows these same expressions with fi ~* 1/4/x. 



lation formatio n redshifts b etween these discretely spaced 
output times C Harker et aL.i200Q ; iGiocoli ct al. 20o3)- The 
open circles and triangles show the corresponding formation 
time distributions from our square-root and constant barrier 
trees. Recall that we do not discretise redshift in our tree, 
so the question of interpolation does not arise. 

For the ellipsoidal collapse model, iGiocoli et all ()2007l ) 
showed that the formation redshift is well-described by 
equation ([7]l if one replaces ujf ~* ^ujf- The smooth 
curves show this with g = 1, 0.707 and 0.55 (short-dashed, 
long-dashed, and dashed), which represent the (constant, 
7 = 0.615, and square-root) barrier predictions. For higher 
values of q, the peaks are located at lower redshifts and 
the widths of the curves decrease. Strictly speaking, equa- 



tion (O only holds for white-noise initi al con ditions. Nev- 
ertheless, as pointed out Lacev & Cofl l)l993l ). it remains a 
reasonable approximation to CDM case. Furthermore, note 
that it provides an excellent description of the formation 
times generated by our trees. However, no choice of q pro- 
vides particularly good agreement with the GIF2 simulation, 
a discrepancy noted by pre v ious authors (.Lin et al. 20031: 
iHiotehs fc Del Popold 120061 : iGiocoh et af] l2007l l. This is 
likely a consequence of the excursion set ass umption that dif- 
feren t step s in the walk are uncorrelated (jSheth fc TormenI 
l2002i ). See |Pan et al] ()2008[ ) and references therein for how 
one might improve on this. 
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Figure 9. Same as Figure [S] but showing only the region around m/M = 1/2. The peak in the simulations (filled symbols) is less 
pronounced than in the merger trees (jagged lines). Open circles show the result of sampling the merger trees at the same redshifts as the 
simulation snapshots: this makes a dramatic difference around 0.49 ^ 0.51, suggesting that the sharp cusp predicted by the theory 
will also be present in simulations with sufficiently closely spaced outputs. The smaller discrepancies further from the peak remain. 



3.3 The mass distribution at formation 

The previous subsection studied halo formation, where for- 
mation was defined as the first time that the mass of one of 
the progenitors exceeds half the total. Therefore, this mass 
can have any value between 1/2 and 1 times the final mass, 
and one can study the distribution of masses at, and just 
prior to, formation. The excursion set constant barrier model 
make s a prediction for this distribution (Nusser & Sheth 
1 19991 ). The mass distribution at formation is expected to 
be 



p(^*)dM= -,/^-^^, where 1/2 «:m<1, (9) 

TT V 2u — 1 U-^ 



and ^ = m/M, and the distribution just before formation is 

where 1/4 ^ ^ 1/2. We have found that, to a very good 
approximation, npip) ^-li^A if o^i^ replaces ^ — > l/4/i 
(solid and dashed curves in Figure [8l respectively). Let ma 
and ruA denote the masses before and after formation, re- 
spectively. Roughly speaking, the symmetry about M/2 in- 
dicates that having a specific ratio of m_B to M/2 before 
formation is equaly likely to having the same ratio of M/2 
to TTiA after formation. 

Although these expressions were derived assuming a 
white-noise power spectrum, they are expected to be rela- 
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Figure 10. The redshift distribution of the most recent major merger, for several final masses. A major merger is defined as one in 
which the minor component has at least 1/3 of the mass of the major component. Filled circles show the GIF2 simulation data, open 
symbols show the corresponding measurements in the same 'snapshot' versions of our trees, and smooth curves show the 'continuous' 
distributions which would be seen with arbitrarily closely spaced output times. 



tively independent of PffcV lSheth fc TormenI (|2004l ') showed 
that they did indeed match numerical simulations well for 
different cosmologies and initial power spectra. Figure [5] 
shows that they also work well for square-root barriers. 

The agreement between the theory curves (smooth 
curves) and our Monte Carlo trees (jagged lines, labeled 
'Cont') is excellent. All panels show that agreement with 
simulation data is also quite good. However, there is a sys- 
tematic discrepancy: the cusp ai ^ — 1/2 appears to be less 
pronounced in the simulation, with corre spondingly lower 
tails. A similar discrepancy was seen bv ISheth fc TormenI 
l|2004h . who suggested that the fact that the simulations 
only provide discrete snapshots in time may be smoothing 
out the peak. By sampling our trees at the simulation snap- 



shots (open symbols, labeled 'Snap'), we have attempted to 
model the magnitude of this effect. Figure |9] illustrates that 
the cusp has indeed been smoothed, but this is a dramatic 
effect only around 0.49 ^ ^5 0.51. The discrepancies fur- 
ther from the peak remain (Figure [8}. 



3.4 The redshift distribution of the last major 
merger 

Mergers of galaxies with similar masses are expected to 
produc e strong short-lived period s of star formation: star- 
bursts ijHernauist fc MihosI Il995h . Recent numerical stud- 
ies suggest that the mass ratio of the galaxies involved 
plays an important role: merger-induced bursts occur 
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when the galaxies have similar masse s llGao et all l2004al : 
ISpringel fc Hernauistll2005l: ICox et al.ll2008l '). Moreover, as 
suggested bv lMaller et al.l (|2006l ). a galaxy's Hubble type is 
strongly correlated with its last major merger. 

Understanding such mergers requires understanding 
the mergers their host haloes undergo. Consider a merger 
(m', m' — m) —> m, with m ' > m — m' . For ease of compar- 
ison with [Parkinson et al.l (|2008h . we will define a 'major' 
merger as one in which (m — m')/m' ^1/3. The filled circles 
in Figure [10] show the redshift distribution of the last major 
merger on to the main branch. (The last major merger does 
not nec essarily happen on the main branch. However, Fig- 
ure 3 of I Parkinson et all ()2008h suggests that, in most cases, 
it does. Presumably, this is because the assembly of haloes 
in recent times is dominated by mergers.) Curves show mea- 
surements in our full trees ('Cont'), and open symbols show 
the result of only sampling the trees at the GIF2 simula- 
tion outputs ('Snap'). For the discretely-sampled data, only 
mergers involving haloes with at least ten particles are con- 
sidered. Note that the anomalously low data point that is 
second from the left in all panels appears to be an effect of 
seeing the tree at discrete snapshots - the smooth curves 
show no such dip. This feature is a l so pre sent in the simu- 
lation analysed by [Parkinson et all (|2008l l: we expect it to 
disappear if more finely spaced snapshots are analysed. 

For high masses, the data from the square-root trees 
peak at about the same redshifts as the simulations; the 
constant barrier, spherical collapse trees peak at lower red- 
shifts. This improvement relative to the spherical collapse 
case is similar to tha t in the modified GALFORM trees of 
I Parkinson" et al.l (|2008l l. However, our square-root trees tend 
to lie above the simulation at low redshifts, and below at 
higher redshifts. The discrepancy with simulation becomes 
increasingly worse at small masses, although it is possible 
that the GIF2 results for our two smallest mass bins are 
not reliable - the high redshift mergers involve haloes with 
few particles. Because we require haloes to have more than 
10 particles, we are likely to miss major mergers once the 
typical mass becomes of this order. 



4 DISCUSSION 

We presented an algorithm for generating merger histories of 
dark matter haloes. This algorithm is based on the excursion 
set approach (Figures (TJ [5] and related discussion) , and can 
handle moving barriers of the sort that are associated with 
ellipsoidal collapse (equation [1} . We illustrated its use by 
generating the forest of trees associated with a square root 
barrier (equation [Sjl. The halo mass function associated with 
this barrier is known to provide a reasonable description 
of halo abundanc es in the GIF2 simulation against which 
we te st our model (iMahmood fc Raieshll2005l :lM oreno et. al.l 
l2008l ). 

Our algorithm produces merger histories which yield 
the progenitor mass functions which are commonly com- 
puted using excursion set theory - demonstrating that our 
approach is nearly self-consistent - over a broad range of 
masses and redshifts (Figures These progenitor mass 

functions are also in reasonable agreement with measure- 
ments in simulations; they are a significant improvement on 
trees based on the constant barrier, spherical collapse model. 



Our algorithm is different from others in the literature 
because it is based on taking discrete steps in mass, whereas 
all others (full N-body simulations included) take discrete 
steps in time. We also used our algorithm to show that while 
the distribution of times at which haloes assemble half their 
mass depends quite strongly on the barrier shape (Figure [7]), 
the distribution of masses at formation does not (Figure [H]). 
Excursion set related formulae for this 'universal' distribu- 
tion (equations [9] and llOp provide an excellent description 
of our merger trees; the agreement with simulations is good, 
but not perfect. 

The algorithm is described in detail in Appendix [X] 
In essence, it requires that one be able to generate one- 
dimensional random walks quickly. Since this reduces to be- 
ing able to generate long strings of Gaussian variables, and 
fast routines for this exist, it is reasonably fast. Significant 
speed-ups can be obtained if one exploits known properties 
of random walks. For example, in the present context, all 
steps in the walk which lie below the current threshold value 
are not interesting (e.g. gray-jagged portions of the random 
walks during the Sm — * S^' jump in Figures [H2|l . If the dis- 
tribution of the number of steps it takes to first exceed the 
current value is known, then one need not generate all these 
steps, one can instead draw a number from this distribution, 
and simply jump this number of steps. Incorporating such 
changes into our algorithm is the subject of ongoing work. 

The excursion-set approach successfully describes many 
properties of the hierarchical growth of dark matter 
haloes. However, with the exception of white-noise ini- 
tial conditions, the progenitor distributions it predicts 
cannot be made consistent with the intuitively attrac- 
tive notion that, in sufficiently small time-steps, merg- 
ers arc binary (Shoth & Pitman 1997; Benson ct al. 2005|; 
iNeistein fc Dekell l2008l : iZhang et all l2008al : |B ensoul |2008| ). 
This accounts for some of the discrepancy between our trees 
and the excursion set based theory curves. A number of 
authors have attempted to alleviate this by relaxing the 
assumption of binary mergers in their merger tree algo- 
rithms. E.g., one of the two algorithms in lSheth fc LemsonI 
l|l999h reproduces the spherical collapse based results ex- 
actly, for arbitrary power spectra. In this algorithm, mergers 
occur between groups of objects rather than just two but, 
as they pointed out , this w a s a rather contriv e d sol ution. 
Somerville fc KolattI l|l999l ), INeistein fc Dekell (|2008l ) and 



Zhang et all (j2008bl l discuss other possibilities. We make no 



attempt to address this issue with our algorithm. 

The appendix also shows that building the main 
branch is straightforward (Figure IA2|I . We expect our 
approach to facilitate studies in volving the mass accre- 



tion history (MAH) of a halo ^Avila- Reese et al 
Nussor & Shcth 1999; van den Bosch 2002; Wcchsler et al 



1998 



20021 : iGao et al.ll200 5: Maulb etsch et al.ll2007l : lstewart et al 



20071 : ICiocoh et al.ll2008a ib). This is the subject of work in 



progress. 

We mentioned in the introduction that halo merger 
histories play an important role in galaxy forma- 
tion models. Models of the merger s of supermas- 



sive black hole s (|Menou et al.l 12001 
I2OO3I : IYoo et al.1 I2OO7I: IVolonteri et al 



induced starbursts (Hernoui st fc Mihosl l 



fiKauffmann_fc Hachnclt 2000l) 
of haloes as their backbone. 



IVolonteri et al.l 
20081 ). merger^ 



995h and quasars 
also have the assembly 
Halo assembly histories 
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also play a key role in studies o f the brightest cluster 
gala xies (iDe Lucia fc Blaizot 2007 ), luminous red galax- 
ies ([Almeida et al.ll2008l: IConrov et al ] |2007l : lMasiedi et alJ 



I2OO7I : IWake et all 



light (jConrov et al 
and the nature of 
last is impor t ant for 



20081) ■ satellites and the intracluste r 



20071 : ISkibba. Sheth fc Martinol l2007i ). 
substructure in galaxy clusters. This 



interpreting the Su nyaev-Zel'dovich 

jHolder et all |2007| ). and strong lensing (|Nataraian et al] 
I2OO7I ) signals^ 

One of the advantages of using Monte Carlo merger 
trees is that one may easily change the underlying cos- 
mology and initial power spectrum. Recently there have 
been attempts to describe spherical coll apse and non-linear 
grow t h in m odificd-gravity theories (jSchafer fc Kovamal 
I2OO8I : llTaszlo fc Bean 2008). In principle, such modifications 
can easily be incorporated into our algorithm. We hope that 
our algorithm will prove useful in some of these studies. 
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APPENDIX A: THE ALGORITHM 

Consider a dark matter halo of mass Mo at redshift Zq 
(Figure rAlf) . Realisations of its merger history may be con- 
structed with our tree algorithm. Below we explain how each 
branch is built and how these branches are connected. 

Al The branch: from random walks to mass 
histories 

The mass growth history of a halo is contained in a ran- 
dom walk (Section l2.3|l . First we discuss the constant-barrier 
(spherical collapse) model and then the square-root barrier 
(ellipsoidal collapse) case. 

A random walk is essentially a collection of steps 
and heights: {so, . . . , Si, . . . } and {ho, . . . ,hi, . . .} (e.g., the 
jagged line in Figure [T]). Consider a horizontal barrier of 
height 5c{z). This line may intersect the walk at several val- 
ues of S. The smallest of such of values, Sm , indicates what 
mass the halo had at redshift z. As z increases, so does 5c - 
and m decreases accordingly, as expected in any hierarchical 
model of halo assembly. 

As we increase the height of the barrier, a subset 
{(So, -ffo), . . . , (Sj, //j), . . . } of the walk is chosen. These 
points contain the mass history (e.g., the dark-filled circles 
in Figure [1} . Every point in this subset has the following 
property: they are higher than all their predecessors. This 
is illustrated by the point at {Sm,5c{z)) Figure [T] it has 
the maximum height in the range Sm ^ S ^ Sm- In other 
words, for any point {si,hi) on the walk to be selected as 
part of the history, it must satisfy the following condition: 

hi > hk, VA; < i. (Al) 

Now we discuss what modifications are necessary when the 
square-root barrier model of ellipsoidal collapse is used (Fig- 
ure [5]). Consider a barrier (equation (3]) with 5- intercept 
y/q5c{z). This curve may intersect the random walk at sev- 
eral places. We pick the smallest of these to find the mass 
at that redshift. As z increases, the (5-intercept increases, 
but the overall shape of the barrier remains unchanged. 
This is how the mass history {{Sq, Ho), ■ ■ ■ , (Sj , Hj), . . .} 
is selected. From this subset we can construct a string 
{yi^Do, . . . , ^/qDi, . . . } of '(5-intercepts', where 

= -^{H, - PVS'S (A2) 

Every point on the history has the following property: its 
5-intercept is higher than that of its predecessors. In other 
words, for a point (si, hi) on the walk to belong to the mass 
history, condition (|A1|I is replaced by 

di> dk, Vfc < i, (A3) 

^ Note that this mapping is ill-defined when 7 > 1/2, because 
these barriers cross. 




Figure Al. A sample merger history tree of a halo with mass 
Mo at Zq. The mass history (black branch) is constructed with 
a random walk. Mass jumps larger than m^^gt are identified, 
and branches are connected there (medium-gray). This process is 
repeated for each branch, and new branches are connected (light- 
gray) until the tree is complete. 

where 

d^ = ^{hi- P^i). (A4) 

For a given cosmology and initial power spectrum one 
may map excursion set variables into physical ones: 
{S{m),5c{z)) ^ {m,z). With this prescription, the mass 
accretion history is obtained: {(Mo, Zo), ■ ■ ■ , {Mj, Zj) . . . }. 
The masses in the history are such that S(Mj) — Sj and 
the redshifts are given by Sc{Zj) = Dj. When barriers are 
constant, Dj = Hj. 

A2 The tree: connecting the branches 

Consider a halo of mass Mo at redshift Zo . We are interested 
in constructing its merger history tree. The first step is to 
draw a random walk with origin at 

so = 5(Mo), ho^ B{S{Mo),5c{Zo)), (A5) 

where B{S,Sc) is given by equation The mass his- 
tory is then collected: {{Mo, Zo), . . . , {Mj, Zj), . . .} (see sec- 
tion IA1[) . In Figure lAll this is represented by the black 
branch. 

The mass history can also be seen as a series of jumps 
in mass: {(Mo ^ Mq, Zo), ■ ■ ■ , {Mj ^ Mj,Zj), . . .}, with 
Mj = Mj+i. Such jumps are interpreted as binary mergers: 
Mj + {Mj — Mj) — > Mj. In practice, we only care about the 
mass jumps above the branching-mass resolution. Denote 
this subset of jumps as {. . . , {mj ^ m'j, z,j), . . . }, where 

mj — m'j > mdust- (A6) 

These jumps are shown in Figure lAll The next step is to 
construct the history of each halo with mass {mj — m'j) that 
falls on to the mass history of the Mo-halo. This is done by 
generating random walks originating from 

So = S{mj ~ m'j), ho = B{S{mj ~ m'j),Sc{zj)), (A7) 
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Figure A2. Main brancli algorithm. The symbols and styles 
are the same as in Figure [2] The dotted vertical line denotes 
S = 5^/2- The main branch is built by following the most mas- 
sive piece at each split; 5^/ > indicates that m' < m/2, 
implying that m — m' > m', whose mass history we must follow. 

Such mass histories correspond to the medium-gray 
branches in Figure ET] The above process is repeated for 
each of these branches: retrieve their mass history, identify 
large enough mass jumps, and attach new branches there 
(light- gray in Figure ET|) . Eventually, all the new branches 
will only involve mass jumps that do not satisfy condition 
(|A6p . At this point the tree is complete and the algorithm 
stops. 

A3 The main branch: a simple modification 

In a merger history tree, the 'main' branch of a halo is ob- 
tained by following the most massive progenitor at each 
mass split. In this section, we show that the mass history 
algorithm described in section IXT1 can be easily modified to 
construct the main branch. This process requires continu- 
ous monitoring of the walk and its associated history, and is 
illustrated in Figure IX2l (compare to Figure [5]). 

Recall that mass decreases with increasing S. For the 
portions in the mass history consisting of jumps where the 
mass loss is less than half, the algorithm is unchanged (e.g., 
the portion with Sm ^ S ^ Sm in Figure fX2)l . Occasionally, 
there are jumps where more than half the mass is lost. Such 
is the case for the Sm — > Sm' jump illustrated in Figure E2l 
with Sm' > Sm/2 (i-e., m' < m/2 and (m — m') > m/2). 
To construct the main branch in that situation, one must 
simply follows (m — m'), not m' . In other words, instead of 
continuing the walk at {Sm' , \/<lSc{z) + P\/Sm'), one must 
continue from {Sm-m' , ^/q^ciz) + Sm-m') (i-e., the dark 
filled circles in Figure 1X2)) . 



© 2008 RAS, MNRAS OOO.ITVITsl 



